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Abstract. The unsteady Micropolar Navier-Stokes Equations (MNSE) are a system of parabolic partial differen- 
tial equations coupling linear velocity and pressure with angular velocity: material particles have both translational 
and rotational degrees of freedom. We propose and analyze a first order semi-implicit fully-discrete scheme for the 
MNSE, which decouples the computation of the linear and angular velocities, is unconditionally stable and delivers 
optimal convergence rates under assumptions analogous to those used for the Navier-Stokes equations. With the 
help of our scheme we explore some qualitative properties of the MNSE related to ferrofluid manipulation and 
pumping. Finally, we propose a second order scheme and show that it is almost unconditionally stable. 



1. The Micropolar Navier-Stokes Equations: Background and Motivations 



The Micropolar Navier-Stokes Equations (MNSE) are a system of time-dependent partial differential equations 
that constitutes a framework to describe the dynamics of continuum media where the material particles have both 
translational and rotational degrees of freedom. Consequently, these equations are very attractive for the dynamic 
description of media subject to distributed couples and polar media in general. 



1.1. The Basic Model. Let us briefly describe the derivation of the MNSE. The mathematical modeling of the 
laws governing the motion of a fluid begins with a description of the conservation of mass, linear and angular 
momentum, which (see [S] or [TJ]) can be written as: 

' ' Dt 

^ (1.1) p^ = diva + pi, 

(1.2) (f + X X u) = pg + px X f + divU + x x diva + CTx , 

where p is the density; u is the linear velocity; cr G M'^^'^ is the Cauchy stress tensor; f is the density of external body 
forces per unit mass; £ is the angular momentum per unit mass; U € M.^^'^ is the moment stress tensor; g represents 
a body source of moments; and (fix )i = eijfctJjfe, where eijfc is the Levi-Civita symbol, i.e., e^j-fc = \{'i—3){j—k){k—i). 
As usual, we denote by D/Dt the material derivative. The physical meaning of the moment stress tensor U is 

y—i analogous to the stress tensor a. In other words, given a plane with normal i^, the vector m = S ■ v is the moment 
^ vector per unit area acting on that plane. 

• '-j Take the cross product of x and ( |1.1[ ) and subtract the result from (1.2 1 to obtain a simplified version of the 

conservation of angular momentum, namely 

^ (1.3) p—=ps + divS + a^. 



Expressions (1.2) and (1.3 1 are usually attributed to Dahler and Scriven (see [1] and [S]) and have been extensively 
used by Eringen (see [8 and _9j) to develop a general theory of continuum media with director fields or, more 
generally, continuum media with microstructure. 

In classical continuum mechanics it is usually assumed that the microconstituents do not possess angular 
momentum and there are no distributed couples. In other words, £ = 0, S = and g = 0. Under these 



assumptions, (1.3) implies that the stress tensor a is symmetric, which is the situation generally considered in 



the literature. These assumptions are appropriate for most practical applications. However, this approach is not 



Date: Submitted to M^AS March 29, 2013. 

2010 Mathematics Subject Classification. 65N12; 65N15; 65N30; 76D99; 76M10. 
Key words and phrases. Micropolar Flows, Ferrofluids, Fluids with Microstructure. 

This work is supported by NSF grants DMS-0807811 and DMS-1109325. AJS is also partially supported by NSF grant DMS- 
1008058 and an AMS-Simons Grant. 



1 



2 



R.H. NOCHETTO, A.J. SALGADO, AND I. TOMAS 



satisfactory (nor even physical) when, for instance, the orientability of the microconstituents plays a major role 
in the physical process of interest. Classical examples are anisotropic fluids, liquid polymers, fluids with rod-like 
particles, ferrofluids, liquid crystals and polarizable media in general. In these cases a precise description of the 
moments and rotations associated to the microconstituents of the material is necessary. 

In the situation described above, the conservation of angular momentum (1.3) needs to be taken explicitly into 
account which, among other things, means that it is necessary to propose constitutive relations for a, i and Z". 
Eringen proposed the following (cf. [71 fT5]): 

^ = Iw, 

where I G M'^^^ is the so-called microinertia density tensor; 

u = {-p + Xdivu)l + ^(Vu + Vu^) + ^,.(Vu - Vu^) + Wx , 

where p is the pressure, I € M.^^^ is the identity tensor, and (wx)ij — £kij^k'i and 

17 = 7o divw I + 7(i(Vw + Vw"^) + 7a (Vw — Vw"^). 

To further simplify the model we will assume that I is isotropic, so that it can be replaced by a scalar j, the so-called 
inertia density. To guarantee that the constitutive relationships do not violate the Clausius-Duhem inequality (see 
[15j). the material constants /x, p,r, 70, 7a and are required to satisfy the following relations: 

3A -t- > 0, > 0, Mr > 0, 7d > 0, 7a + 7<i > 0, 

(1.4) 

870 -t- 27d > 0, -(7a + 7d) <ld~la< {la + Id)- 

As a final simplification, we will assume that the fluid is incompressible and has constant density. 

Let C M'' with d = 2 or 3 be the domain occupied by the fluid. Replacing these constitutive relationships 
into and pTsI , we arrive at the MNSE, 

Ut — (i^ + ^'r)Au H- (u • V)u + Vp ~ 2i/rCuiiw + f , 
(1.5) { divu = 0, 

jWf ~ (ca + Crf)Aw -I- j(u • V)w — (co + Q — Ca)'Vdivw + AUrV/ — 2i^rCurlu + g, 

where we implicitly redefined the pressure as p^^p, and the constants ly, Vr, Ca, Cd and cq are the kinematic 
viscosities (i.e. /i, ^r, 7a , Id and 70 divided by p, respectively). This system is supplemented with the following 
initial and boundary conditions 

, , u|t=o = uo, w|t^o = Wq, 

(1.6) 

u|anx(o,T) = 0, w|as2x(o,T) = 0. 

The reader is referred to [TS] for questions regarding existence, uniqueness and regularity of solutions to ( 1.5 1-( 1.6 ) 
and related models. The purpose of our work is to propose and analyze numerical techniques for this problem. To 
simplify notation, in what follows we will set 

(1.7) V„=V + Ur, Ci=Ca+Cd, C2 = Co + Cd - Ca, 

and we will assume that ci, C2 > (see for instance |15j ) which is consistent with the thermodynamical constraints 
0. 

The MNSE can be regarded as a building block of models that describe the physics of polarizable media. For 
instance, Rosensweig (see [55]) described the behavior of ferrofluids subject to a magnetizing field h with the 
MNSE and 

{f = Mo(m-V)h, g = Momxh, 
1 
mt — aAm + (u • V)m = w X m — — (m — >foh) mil, 

where m denotes the magnetization field and ,T > 0, a > 0, xtq > are material constants. The magnetizing field 
is assumed to obey the Maxwell equations. The reader is referred to [T] for an analysis of this model. 

In addition to applications in smart fluids and polarizable media, there has been a growing interest on the 
MNSE in other areas. For instance, they have been used to describe collisional granular flows, where the size of 
the microconstituents is comparable to the macroscopic scale ([H]) and the frictional interaction between particles 
is not properly modeled by the classical equations of hydrodynamics. Another application is the modeling of micro 
and nano flows ([20]). where again the size of the microconstituents is comparable to the "macroscopic" scale and 
the rotational effects cannot be neglected. 
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Figure 1. Idealized configuration of a ferrofluid pumping experiment. A planar duct with a 
solenoid that generates a uniform magnetizing field h = hoi. Since, g — Hom x h (see (1.8)), 
it will produce torque in the regions where h and m are not collinear. In a real ferrofluid the 
magnetization vector field m would evolve through the channel satisfying the evolution equation 
( 1.8 1 and will try to align with the magnetizing field. However, and as part of an idealized setting, 
we will assume that the magnetization profile m depends only on the j/-direction. 



The key points of this paper are organized as follows. Section 1.2 introduces a very simple experiment (ferrofluid 
pumping) as a motivation for the analysis and numerical implementation of the MNSE. In Section 2.1 we recall 
the basic energy estimates and existence theory for the MNSE. Paragraphs |2.2| and |2.3| introduce the notation 
and the basic tools required for the analysis of the numerical scheme proposed later in Section [3j Error estimates 
for the linear and angular velocities are derived in Section [4. 1[ and error estimates for the pressure are derived in 
Section [4.3[ We present a formally second order scheme in Section [5j and show that it is almost unconditionally 
stable, i.e. it is stable provided the time step is smaller than a constant dependent on the material parameters, 
but not on the space discretization; see (5.51 for details. Finally, in Section |6j we provide numerical validation of 
the error estimates derived earlier. 



1.2. Potential Application: Ferrofluid Pumping by Magnetic Induction. To illustrate the differences 
between the MNSE and the classical Navier-Stokes equations here we propose a setting by means of which it is 
possible, at least theoretically, to generate fluid motion via a well designed forcing term in the equation of angular 
momentum. This example is inspired by j25i . where a ferrofluid is pumped by the actuation of a spatially-uniform 
sinusoidally time-varying magnetizing field. Another pumping strategy, this time based on a magnetizing field 
that is varying in space and time, is proposed in [TO . 

The idealized setting that we shall consider is depicted in Figure [l] We assume that our domain is a planar 
duct of unit height and length L > 1, which is wrapped by a solenoid that generates a uniform magnetizing field 



h = hoi, where hp is just a positive constant. From ( 1.8 1 we infer that f = 0, since the magnetizing field is constant 



in space. As part of our idealized setting, we disregard the evolution equation in (1.8) for the magnetization field, 
and set m to be constant in time and depend only on the vertical variable y, i.e., 

m = mo (cos di + sin 9j) , 



where toq is just a positive constant, and 9 = 6{y). Using (1.8) we get 



(1.9) g = -/ior7io/iosin6'(y)K. 

As reference configuration we will consider a linear interpolation between the points (0,7r/2) and (1,— 7r/2), 
that is 
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Figure 2. Plot of the function 9i{y) (dotted line), and the family of functions {di{y)}'l^2 (solid 
lines). These are used to induce a force in the angular momentum equation. The function 9i is a 
linear interpolation between ±7r/2 and 9i, for i = 2, . . . ,7 are small perturbations of it. 



g2(y) 



9i(y) 
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Figure 3. Velocity profiles obtained with the forcing terms {gi}I=2 (solid lines). For comparison 
the velocity profile obtained by using gi is also shown (dotted line). The figures show that it is 
possible to generate linear velocity via appropriate actuation in the angular momentum equation. 
Notice that, although it is not dramatically different from the others, the forcing term g7 induces 
motion in the opposite direction. 



As perturbations from this reference case we consider, for z = 2, . . . , 7, 

7r(480a;5 - 1200a;4 - Ax^{i^ + IQi - 275) + Qx^{i^ + IQi - 75) - - 5(2i - 7)) 
'^y' " 2(*2 + 10* - 35) ■ 

A plot of these functions is provided in Figure [2] Notice that they aU satisfy 6li(0) = 7r/2, 61^(1/2) = and 
9i{l) = — 7r/2 which we require to model a magnetization field that is perfectly aligned with the magnetizing field 
at the center of the channel, but is perpendicular to it at the top and bottom walls. 

We assume the fluid is initially at rest, the boundary conditions for the upper and lower part of the duct are no 
slip, and for the left and right sides of the duct we consider open boundary conditions. We apply the magnetizing 
field linearly in time, that is we set h = ho{t/T)i. We let L = 1, and the material constants he i' = i^r — ^, 
Ca ^ Cd = cq = 1, and j = 1. We use a Taylor-Hood finite element discretization of 40 elements in the horizontal 
and vertical directions, and a time-step r = 1/50. The numerical scheme used for this example is the first order 
method discussed and analyzed in this work. Figure [3] shows the velocity profiles at time t = T and x = 1 obtained 
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by setting g as in (1.9). These results are stable (in the sense that they do not change) with respect to the spatial 
and temporal discretizations, and length of the channel. However, as it would happen with any physical model, 
these results can be sensitive to changes in the constitutive parameters. A discussion about the possible influence 
of the constitutive parameters on the pumping phenomena goes beyond the scope of this paper (see for instance 

The results in Figure [3] give an idea about the kind of forces that are necessary in a real ferrohydrodynamic 
setting, in particular in the case of a spatially uniform and sinusoidal in time magnetizing field as in |25) . The main 
observation here is that small variations of the forcing term can yield quite different flow regimes, including flow in 
the opposite direction, this feature is observed in experiments (cf. [26] ). Finally, the reader should be reminded 
that this is just an idealized setting which illustrates the main pumping mechanism. In real ferrohydrodynamics 
we cannot set the value of magnetization m as we please because m is actually determined by the evolution law 
in ((L8|. 



2. Notation and Preliminaries 



We shall consider system (1.5 1 in an open, bounded, simply connected domain f2 C R'^ with d = 2, 3, with a 
smooth boundary dfl, for a finite interval of time (0, T), and we will denote = fl x {0,T). We use the standard 
Sobolev spaces W^(r2) for < fc < oo and 1 < g < oo that consist of functions / G L'^{n) whose distributional 
derivatives of order up to k are also in L'^{^l). To simplify notation, we set H''{il) = (i7), and denote the 
closure of C^(ri) in if^(r2) by HQ^fl). We denote with bold characters vector valued functions and their spaces. 
The scalar product in L^(ri) and L^(17) are indistinctly denoted by (•, •). The subspace of functions in L'^{n) with 
zero mean is denoted by LQ{n). Whenever i? is a normed space, we denote by |j • \\e its norm. The space of 
functions </> : [0, T] -J> S such that the map (0,r) 9 t ^ ||(/)(i)llE e R is L'-integrable is denoted by L'i{E). 

We shall make repeated use of the following integration by parts formula for the curl operator: 

(2.1) (curiw,u) (w,curJu) Vu,weHj(f]). 
In addition, we recall that the following orthogonal decomposition of HQ(r2) 

||Vu||2, = \\curlu\\l2 + ||ciivu||2,, Vu e HJ(1]) 
holds true (provided D, is bounded and simply connected, see for instance [11 ) which implies 

(2.2) \\cmlu\\l2 < ||Vu||2, Vu e UHn). 

We use the following two classical spaces of divergence- free functions (see for instance [53]) 

H = {u e L^{fl) I divu = in f7 and u • = on dfl} , V = Hj(r2) n H. 

Henceforth C denotes a generic constant, whose value might change at each occurrence. This constant might 
depend on the data of our problem and, when discussing discretization, its exact solution, but it does not depend 
on the discretization parameters or the numerical solution. We denote by Cp the best constant in the Poincare 
inequality, i.e., 

||u||l2 < Cp||Vu||l2 Vu e Ul{n), Cp « diam{n). 
We will use, as it has become customary, the following trilinear form 



6(u,v,w)=V'/ u'v^.w^dx, u, V, w e Hj(ri), 



which, as it is well known (cf. |^ ) , is skew-symmetric whenever the first argument belongs to V. In addition, we 
shall use the following, also well known, inequalities (see [TT]): 

(2.3) 6(u,v,w) < q|Vu||L2||Vv||L2||Vw||L2, Vu, v, w e Hi(l]), 

(2.4) &(u,v,w) < C||u||l~||Vv||l2||w||l2, Vu e H^{n),\f^r e Hj(17),Vw e L^il), 

(2.5) 6(u,v,w) < C||u||l2||Vv||l2||w||l-, Vu e L^{n),yv G Hj(fi),Vw e U^{n), 

(2.6) 6(u,v,w) < C||u||l2|1v||h2||Vw||l2, Vu e L^n),y^ e H2(f]),Vw e Hj(fi). 
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2.1. Energy Estimates and Existence Theorems. The stabihty and error analysis of the scheme that wiU be 
proposed in Section [3] is based on energy arguments. Therefore, to gain intuition, let us briefly describe the basic 
formal energy estimates that can be obtained from (1.5). Multiply the linear momentum equation by u and the 
angular momentum equation by w and integrate in J7. Adding both ensuing equations, we obtain 



1^ 
2dt 



lL2 



2 2 2 2 / \/\/\ 

Vu||l2 + ci ||Vw||l2 + C2 ||djVw||L2 + 4i^r ||w||l2 = 4iyr(curiu, w) + (f, u) + (g, w) 



where the parameters Vg, c\ and ci were defined in (1.7 1. Repeated applications of Young's and Poincare's 
inequalities yield, after integration in time. 



|u(t)|| 



L2 



(2.7) 



-Jl|w(t)||^, 



■ [ ||Vu(s)||l2 ds + ci I ||Vw(s)||l2 ds + 2c2 j ||divw||L2ds 
Jo Jo Jo 



< 



< c: 







'"f(^)lli' ' 



2 

L2 



ds 



l|2 

|uo|Il2 



jIIwo||l2 yt<T, 



V Cl 

This formal energy estimate suggests that solutions to (1.51 are such that 
(2.8) ue L°°(H)nL2(v)^ weL°°(L2(17))nL2(Hj(l])). 

To obtain an estimate on the pressure, we use a well-known estimate on the right inverse of the divergence 
operator (cf. [HI [6]), i.e., 

[q, divv) 



(2.9) 



l3\\q\\L^ < sup 



I Hi 



From (2.9 1 and the linear momentum equation in (1.5) we get 

i-T 



hds< 



{\\nt{s)\ 



2 

L2 



|Vu(s) 



2 

L2 



|Vu(s) 



L2 



|Vw(s)||2, + ||f||2,)ds, 



SO that, to obtain an estimate on the pressure, we must assume u G L*(Hj(r2)) and, in addition, we need an 
estimate on the time derivative of the linear velocity at least in i^(L^(i7)). This is standard for the Navier-Stokes 
equations. To obtain it we differentiate with respect to time the equations of conservation of linear and angular 
momentum. Repeating the steps used to obtain (2.7 1 we arrive at the desired estimate. 
The existence of weak solutions can be summarized as follows. 

Theorem 2.1 (Existence of weak solutions). Let f , g G iy^(L^(57)), Uq G H and wq G L^(il). Then there exist 
(u, w,p) G i°°(H) X L°°(L^(f2)) X 'D'{flT) satisfying ( 1.5 ) in the sense of distributions. Moreover, u and w satisfy 
the energy estimate (2.7). 



Proof, see [T5l Theorem 1.6.1]. 

Just like for the Navier-Stokes equations, uniqueness of solutions of the MNSE is an open issue. 



□ 




2.2. Time Discretization. We introduce A' > to denote the number of steps, define the time-step as t = 
T/K > and set t'' ^ kr ior 0<k < K. For cj) : [0, T] E, with E being a Banach space, we set = (/((i*^). A 
sequence will be denoted by (j)'^ = {0'°}j._q and we introduce the following norms: 

||0ni£~(iJ)=„max^ll0'lU, Ul\iHE) = 

We define the backward difference operator 

(2.10) 5/ = c/)*^ - </>'=^\ 

and set 5^(1)'' = 5(5/) = cf)'' - 20'^- ^ + (j)^-^. 

Finally, recall the following discrete Gronwall inequality. 

Lemma 2.1 (Discrete Gronwall). Let a7 , b'^ , and 7"^ be sequences of nonnegative numbers such that rjk < 1 
for all k, and let go > be so that the following inequality holds: 

K K K 

+ T 6fe < r ^ 7fcafc + r ^ Cfc + 50- 
fc=o k=o 



/ , Ok 

k=0 
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Then 

Oif + r ^ 6fc < r ^ Cfc + go exp r ^ (ikjk , 
fc=o V fc=o / V fc=o / 

where ak = (l - r^k)'^ ■ 

Proof. See mill]. □ 



2.3. Space Discretization. To construct an approximation of the solution to (1.5 1 via Galerkin techniques 
we introduce two famihes of finite dimensional spaces, \yh\h>o and {Q/i}h>o with V/i C Hj!j(ri) and, Q/j C 
H^{Vl)C\L'^{Vl). The space will be used to approximate the linear and angular velocities and Q/i to approximate 
the pressure. We require that these spaces are compatible, in the sense that they satisfy the LBB condition 

/r, 11^ • r {qh,divvh) ^ 

(2.11) mf sup - — - — — - — - — > p , 

0#9h SQh o^-Uh 6 V, J 1 9/1 1 1 L2 1 1 ""^/i 1 1 L2 

with /?* independent of the discretization parameter h. In addition, we require that the spaces have suitable 
approximation properties, in other words, there exists a, I € N such that for m G [0, 1], 

inf \\q - qhh- < Ch'^\\q\\H^. Vq E F"(0) D Ll{n), 

inf (||v-i;,.||l2 +/i||v-t;^||Hi) < C/i™+1||v||h™+i , Vv e H'"+i n Hi(r!). 
Lastly, we assume that the velocity space Wh satisfies the following inverse inequality: 

(2.12) ||u;,||l- < C7^(/i)||u,,||hi Vu,, e V,,, 

where = (1 + | log(ft,)|)5 if d = 2 and ^{h) = h^^ if d = 3. References [101 E] provide a comprehensive list 
of suitable choices for these spaces. 

For a.e. t G [0,T] we define the Stokes projection of {\i{t) , p{t)) as the pair {\ih{t),ph{t)) G V/i x Q/j that solves 

J {VvLh, Vvh) - {ph, divvh) = (Vu, \/vh) - {p, divvh) Mvu G V/i 
I {qh, divxih) = {qh, divu) ^qh G Qh ■ 

In addition, we define the elliptic- like projection of w(i) as the function Wh{t) G that solves 

Ci(Vw/i, Vw/i) + C2{divwh, divujh) + 4:i'r{^h,^h) = 

= ci(Vw, VLOh) + C2(divw, divujh) + 4;/,-(w, u;,i) Vu;/i G Wh ■ 

The properties of the Stokes and elliptic- like projections are summarized in the following; see for instance |12| . 

Lemma 2.2 (Properties of projectors). // (u,w,p) G [L°°{U^{fl) n Hi(f]))]2 x L°°{H^{n) n Ll{n)), then the 
Stokes and elliptic-like projectors are stable in dimension d < 3, i.e., 

(2.13) |lUft||L-(L=°nWi) + l|wh||L-(L°°nWi) + \\Ph\\L-^{H^) < C (||u|1loc,(h2) + |b||L~(Hi) + || w|l (h^)) . 

//, in addition, (u,w,p) G [L°°(H'+i(fJ)nHj(n))]2 x L°°(i7'(f7)nL^(f2)), then the projections satisfy the following 
approximation properties: 

||u-U,,||ioc(L2) -^/l||u-U,i||^oo(Hi) +^||P-P/i|1l-(L2) < C/l' + ^77(u,p) 

(2.14) _ 

||w - Wft||i^(L2) + h\\w - w,j||ioo(Hi) < C/i'+^$(w), 

where 77(u,p) = |1u|1^„„(ht+i) + |b|l^oc,(jjr) and ^{w) = ||w||^^(jjt+i) . 

We introduce the trilincar form bh : [HQ(f2)]'' — > M, 

bh{u, V, w) = b{u, V, w) + -(divu, v • w), V u, v, w G Hj(0), 

and recall that it is consistent, i.e., 6/i(u, v,w) = 6(u,v,w) whenever u G V, and skew-symmetric 

6/i(u, V, v) = 0, 
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for all u, V G HQ(fl). This form satisfies estimates similar to (2.3)-(2.6), namely 

&/i(Uh, V,,, W,j) < C|lVUh.|lL2||Vv/jL2|l'^W,J|L2, Vu,,,V,j,W/, e V;,, 

^'^'^^^ &/i(u/i,v,w,0 < C||u,jL2||v||H2||Vwft||L2, Vuft,w,, eV,„ VveH2(17), 

and 



(2.16) bh{uh,Vh,Wh) < C||uh||L2||v/i||(LO°nwi)l|Vw,,||L2, Vuh,v,i,w,i e Yh- 



Since, by assumption, the space V/j satisfies the inverse inequality (2.121, then for d = 3, 



bh{uh,^hi^h) <Ch H|u/iIIl2|1Vv,,||l2||Vw,i||l2 Vu,j, v,j, w,^ e V?, 

(2.17) 

^/i(uh, v/j, w/i) < C/i HI'^U'!||L2||Vv/i||L2||wft.||L2 Vu/i, v,i, e V?i 

3. Description of the First Order Scheme 

To the best of our knowledge, the only work that is concerned with the construction and analysis of a scheme 
for the MNSE is [19], where a fully discrete penalty projection method for this system is developed and analyzed, 
and a suboptimal convergence rate is derived. Our scheme instead possesses optimal approximation properties 
and requires the solution of a saddle point problem at each time step, which can be done efhciently. However, it 
can be easily modified to decouple the linear velocity and pressure via an incremental projection method, while 
maintaining optimal orders of convergence. For brevity this will not be included. 

Let us now describe the scheme. The scheme computes {U^,W^,P,[} C x meant to approximate, at 
each time step, the linear and angular velocities and the pressure. Wc initialize the scheme by setting 

(3.1) (U°,P°) = (uO,p°), W°=w°, 

that is, we compute the Stokes and elliptic-like projections of the initial data. 



Remark 3.1 (Initialization). The initialization step (3.1) requires that the initial data is regular enough so that 



the projections are well defined, which from now on we will assume. If this is not the case, (3.1 1 must be modified 
and, say, take L^-projections. The analysis below must be accordingly adjusted to take this into accoimt (cf. |14)). 

After initialization, for fc = 1, . . . , i^T, we march in time in two steps: 
Linear Momentum: Compute (U^,P^) e V/i x Qh, solution of 

(3.2a) +iy,{VUlVvh) + b4ut\VliJh) - {PtdivvH) = 2j,r{curlW'^-\vh) + (i^Vh) , 

(3.2b) (g^,divU^) =0, 

for aU Vh e Wh, qh e Qh- 

Angular Momentum: Find G Yh that solves 
(3.3) J +ci(VWtVc^,)+j5,(Utwta;,0 + 



for all ujh G Yh- 



C2{divWldivuJh) + Ai'ri^ti^h) = 2i^^(curJUtu;,,) + (g^a;^) , 



Notice that we have decoupled the linear and angular momentum equations by time-lagging of the variables. 
This scheme is unconditionally stable, as the following result shows. 

Proposition 3.1 (Unconditional stability of the first order scheme). The sequence {U^, W^,P^} C [V/J^ x Q/j, 



solution of (3.2 )-( |3.3[ ), satisfies 

K K 



I|U^IIl2 + (j + 4^rr)||W^||^, + (I|5U^IIl2 + J\ml\\h) + E ^ (H|VU;;||^, + rcill VW^II^,) 
(3.4) ^ ""^^ 

+ 2j2rc,\\divWl\\l. <j2r{^\r\\h + ^l|g''1l^2) + ||U?J|^. + (j + A.rr)\\Wl\\h . 

k=l k=l \ ^ J 
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Proof. Set = 2rU^ in ( |3.2[ ) and iv^ = 2rW^ in (3.31, respectively, and add the results. Use the identity 
2a{a — b) = — + {a — bf, the integration by parts formula (2.1 1, estimate (2.2) and Young's inequality to 
obtain 

IIU^IIl^ + (j + 4i^rr)|lW^||2, + ||5Ufj|2, +ji|8Wfj|2, +rHlVU!;||L +rci||VW^||2, +2rc2||divW;;||^, 

(3.5) 



<\\Vt'\\i^ + ij + ^'^rr)\\W^-'\\i 



f fc||2 



L2 



Cl 



Adding over k we obtain the desired estimate (3.4 1 



□ 



4. A Priori Error Analysis 



Here we perform an error analysis of scheme (3.2|-(3.3| and show that this method has optimal convergence 



properties. The analysis is based on energy arguments and hinges on the unconditional stability result of Proposi- 
tion [STT] The arguments used are rather standard for the Navier-Stokes equations, the main novelty and difficulty 
being the coupling with the angular momentum equation, which requires lengthy and careful computations. 



We shall assume, for the sake of simplicity, that the solution to ( 1.5 1-( |1.6[ ) satisfies: 

(4.1) u,weCi([0,T],H'+i(O)) and Uu,wu e L\[0,T],L^n)). 

These assumptions will be enough to derive optimal convergence rates for the linear and angular velocities. If we 
want to do the same with the pressure we will require the additional regularity: 

(4.2) uu,w« eC([o,r],H'+i(f^)). 

These assumptions are standard in the error analysis of incompressible flows (cf. [IT]). 

The first step in the error analysis is to analyze the consistency of the method. To do so, we proceed as it is 
customary in the analysis of evolutionary problems (cf. |24j ) and split the errors 



P 



into the so-called interpolation and approximation errors via the Stokes and elliptic projections of ^2.3 
(4.3) 



I.e. 



Ph ■ 



Tjfc 



ck 



pI 



pk 



The interpolation errors (S"^, R'^, r'^) are controlled by means of (2.14), so that the next step is to derive an energy 



estimate for the approximation errors (E^, e^) which is a slight variation of that one obtained for (U^, W^, PJ^) 
in dslf. 



4.1. Error estimates for the Linear and Angular Velocities. The approximation errors (E^,f^,ep satisfy 
the following energy identity: 



k\\2 



7^k-l\\2 



-fe-l||2 



(4.4) 



i.,+2Tv,\\yK\\h 

6 



2TCi||Vf,t||^2 + I|5eS||2 



2TCy\\div£, 



k\\2 



with 



^3 
As 

A, 



2r 6„(U^i, E;;) - 2r 6,(u^ u^ E^) 
2jT bn{\5l Wl ~ 2jT 6,(u^ w^ 4^) 
4Ti/,.(curJw''' - curiW^"\ E^) 
4T:/^(curiE^£:;5:) 
-2(5S^EfJ-2J(5R^£,t) 
2r(7^^;,EfO+2JT«,4% 



where TZ'^ and TZ'^ arc integral representations of Taylor remainders (see for instance [H]), i.e. 



(4.5) 



K = - t it"-' - s)utt{s)ds and K ^ - t {t"-' 

T Jtk-l T Jfk-l 



,s)wtt(s) ds. 
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The main difficulty, and our focus from now on, is to estimate the residual terms Ai, i = 1, . . . ,6. 



Theorem 4.1 (Error estimate on velocities). Assume (4.1), then 
(4.6) 

whenever 



(4.7) 

where M and A4 satisfy 
(4.8) 



EIIloc(l^) + ||f1lL»(L^) + h (||VEnU2(L.) + \\V£^\\l2^l^)) < C{t + h'+^) 



r < — with K ■ 
K 



vo' ci ' ci ' ci ' 



sup(||Vuft|lL3 + |1u,,|1l~) + sup |u|2 < M < oo , 
sup(||Vw/i||l3 + ||w,j||l~)^ + sup |w|^ < M < oo . 



Proof. It suffices to provide bounds for the terms Ai above and employ the discrete Gronwall lemma. To begin 
with, notice that 



(4.9) 
and 
(4.10) 



1 ,,k 



,k — l T7\k 



- bh{El<,^h) ~ 5,(S^w,';,u;,) - bh{ul£t,^h) Va;, G V, . 
Set Vh — in ( |4.9[ ). Since bh is skew-symmetric the last two terms vanish, and we can rewrite Ai as: 

Ai = -2t ^5u^ u^ Ej;) - 2r 6„(u'^-\ S^ Ef,) - 2r 6^(S"-\ E^;) - 2r 6;,(E^-i, E^) 

The functions hu^ and u'^~^ are solenoidal so that the consistency of bh yields control on An and A12: 



An = 2r6,(5u^Etu'^) < 2t||5u^-||lHI VE^||l. ||u1|l- < ^IIVE^H^, + ^HSu'^^H^, 



= -2r6,(u^-^S^E^) < 2t IIu'^-iL^ || VE^||l. ||S1|l^ < ^IIVE^H^, + — HS'^H^,, 



where we have used (2.5) and ( |2.4[ ). By (|4.1|), we deduce 
(4.11) 



l|5u' 



< T 



The terms Ai^ and A14 can be estimated via (2.16) as follows 

9Mt2 



Ai3 + Ai4<^||VE;;||2 



L2 



Sfe-i||2^^9Mr2||E^i„2 



L2- 



Set = 2t£^ in (4.10). We rewrite A2 as 

^2 = -2Jr^u^R^£;J) -2Tjfo,(Etwt£,t) -2rJ6;,(S^w;;,£,^) - A21 + + A23 



Since u is solenoidal the bound on A21 proceeds as that of A12, whereas (2.16) gives control on A22 and A23: 



(IIR'^IIL 



|E,';||2. + ||s^-"2 



7 - ci 

The bound on A-^ begins by noticing that w'^' — Wj^^^^ = Sw*' + R*"'^^ + ^^^^ ■ The integration by parts formula 



(2.1) then yields 



A3 = 4r^y^(5w^ cur7E^) + 4ri/.r(R''"\ curiE^;) + 4Ti^^(£:^'"\ curiE^) 
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whence 



^3 < 4ri.,||5w'=||L.||VE^||L2 +4Ti.,||R'^-i||L2||VE^|| +4Ti.,||£;J-i||L2||VE^||L2 



<^I|ve;;||L + ^I 



L2 



7fe-l||2 



The term ||5w''|jL2 can be bounded similarly to (4.111. 
The bound on follows the same lines as those of A3: 

A4 = ATVr{^\ curl Et)+ATVr{K, curl £f,) < 

The last two terms A5 and Aa can be easily bounded as follows 

J 

and 

J/„T. .... 9C?.T 



fell 2 
L2 



Cl 



A. = -2(5S^E;;) - 2J(8R^£,^) < ^||VEfj|2, 



P||§Sfe,|2 



A6 = 2r(7^^,E^) 



2rX7^t,f^■)<i-|lVE^IlL^ 



p 



l7^ 



L2 



fe||2 

uIIl2 



ClT 

7 ' 

ClT 

7 ' 



28i^,2r 

Cl 



I^||§Rfc„2 
ClT 



L2 ) 



(4.12) 



The interpolation errors are bounded by (2.14) which, in conjunction with (4.1), also implies 

IISS'^IIl^ +ft||5VS'=|lL2 <CTh'+^ij{ut,pt) 
||5R'''||l2 +/i|15VR'=||l2 < Ct h}+^£,{MVt). 



Assumption (4.1) also gives an estimate on the truncation errors TZ^^ and TZ: 



k 

w 1 



(4.13) 



\'<\\h < 



\uu\\hdt , IKII^. < 



3 7t 



l^ttWh dt. 



Inserting the estimates above for Ai, 1 < i < 6, into (4.4), summing in k and application of Gronwall inequality 
concludes the proof. □ 



Remark 4.1 (Smallness assumption on r). Condition (4.7) does not depend on the space discretization parameter 
h. It does depend, however, on the constants M and A4 defined in (4.8); this is standard for Navier-Stokes. In 
addition, this estimate depends on the quotients v^jv^^ and v'^/ci, which gives an indication of how strong the 
coupling between linear and angular momentum is. 

4.2. Error Estimates for the Discrete Time Derivative. When dealing with the Navier-Stokes equations, it 
is well-known (see, for instance, [T^]) that in order to derive optimal error estimates for the pressure in £'^{L?'{Vl)) 
one must first obtain estimates on the discrete time derivative of the velocity, which is the main reason for the 
additional regularity requested in ( |4.2[ ). Our analysis is no exception, and this is additionally complicated by the 
fact that we must obtain error estimates for the derivatives of the linear and angular velocities. However, it is 
important to point out that it is possible derive an error estimate. 

Applying the increment operator 5, defined in (2.10), to the equations that govern the approximation errors and 
proceeding as in the proof of Proposition 3.1 we conclude that the discrete time derivatives t^^SEJJ and r^^hSJ^ 



satisfy an energy identity similar to ( |4.4| , namely, 

||r-i5Ef;||^, + (^^ + 8^.r)||r-i54^||^. - ||r-i5E^i||2, - j||r-i5£:ri||2, + 2^„r||r-i5VE^||2, 
(4.14) 5 



+ 2cir||r-i5Vf,^f + ||r-i5XllL^+j|k"'52f/tllL+2c2r||r-Miv5f^-||^, 

1=1 



where 



Fi - 2 6;,(U^\ Ut t-iSE;;) - 2 fe, (u^ u^ T-i5E^) - 2 6„(U^2, U^i, t-iSE;; 

Fy, = Avr{curlhvf'' ~ curi5W^-\T-i5Ef;) +4i/^(curJ5u'= - curiSU^ r-^Sf^) 
F4 = -2r-l(52S^r-l5E^) -2jT-i(52R'=,T-i5f^) 
F5 = 2(57ei;,r-i5E;;) +2j(5<,T-i5f;J). 



2fe,(u«-\u«-\r-i5E;;: 



,fe-i „fe-i 

,/c-l ,„fe-l 



2j6Ju«-\w«-\T-i5£-,t) 
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A bound on these terms then yields a bound on the discrete time derivatives. This is the content of the fohowing 
result. 



Theorem 4.2 (Error estimate for the discrete time derivatives). Assume (4.2). // 



(4.15) 

are sufficiently small, then 
(4.16) 



/i-i/2||e;;||,o.(l.) and h-'/^£l\ 



5EI 


+ 




T 




T 



< C{T + h}). 



^(L2) 



Proof. In analogy to Theorem 4.1 it suffices to bound the residual terms The proof is rather technical and 

tedious, and consists of careful manipulations of these five terms. Take the difference of (4.9) for two consecutive 
time-steps, which allows us to write Fi as the sum of six terms {Fu}^^^: 



Fii = 


-2 6(5u^ 


u^T-l5E|;) + 


2b{5u''-\u''-\T- 


-iSEfj 


Fl2 — 


-2fo(u'=-^ 


,S^r-l5E^)- 




-i5E^) 




-2fo^(S'^- 




+ 2bh{S''-^,uf-\ 




Fl4^ 






+ 2 6^(E|;-2,u|;-i 




Fl5 = 






-2 6,(Er^E^i,r-i5E^) 


-^16 = 




-i,Etr-i5E^ 


+ 2 6,(u^2,E^i 


,r-i5Efj 



Using the linearity and skew-symmetry of the trilinear form, these six terms can be appropriately rewritten and 
bounded using ([23]) -([ZS]) and ( |2.15| -( |2.17| to get 

2 6(5u^ T-'ml, Su*--) + 2 6(52u^ T-^SEt u'^-^) 



11 



< ^||r-i5VE^||^. + — IISVu'^ll^. + '^\\r-'6WE' 
14 i/„r 14 



,fe — 1 _ — ISTJife 



^1 



Fi2 = 2 biu^-^T-'dEl, 58") + 2 6(5u'^-\ T-^SE*, S 



- 14 II "IIL II IIL 



14 ' 



r-i5VEf;||^. 



c 



L2 



ife-l|i2 



3fc-2 



^^13 = -2 6„(5S'=-\utT-i5EfJ -2&„(S 

M 

iirts^^-^ii;., -I- 

14 

^^14 = -2r6^T-i5Ef-\utr-i5E;;) - 2 6;,(Ef^ 

;ife-l||2 



<^||r-i5VE^||^. + -^||5S'=-i||2, 

ri/c-l 



--i5VE'=ii2 



fe||2 



|5V<-i||^.| 



/c-2 



5<,r 



-i5Ef;) 



< — IIT 

- 14 II 



r-'?>El- ,|L2 



14 



|r-i5VE^||2, 



C 



|8Vu;;||^,| 



VE 



fe-2||2 



L2 I 



^^15 = 2r6,(r-i5Er\Er\r-i5E!;) < ^^^^r (||r-i5VEri||^ 



Fi6 = -2b,i8nl-\Elr-'mt) < ^||r-i5VE^||^, 



^1 



L2 ^ 

L2||VE^||L2 



^-i5VEfeii2 



5yufe-l||2 ||y7T7fc||2 



Similarly, applying 5 to (4.10), F2 can be expressed as the sum of five terms {^2j}j=i- 

i^2i = -2J6„(u^R^T-l5f^•) +2j6„(u^-i,R'^-i,r-i5£:;J) , 
F22 = 2j6,(Et£^T-i5f,t) - 2J64E^^f^^r-l5£:;^) , 
i^23 = -2j6;,(Etwtr-i5£^^) +2j&,(E^\w^\r-i54^) , 
i^24 - -2J6;,(S^wtr-l5£^^) +2j&,(S'=-i,w|;-\r-i5£^^) , 
F25 = -2j6;,(ut£:^r-i5£:;J) +2j6;,(u^i,£:^i,T-i54^) . 
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We now bound each of these terms separately 

F21 - -2J&„(5u^R^T-l5£^) +2J64u'=-^5R^T-l5£,t 



12 ClT 



5Vu'=||2,||VR'=||2, 



ClT 



12 ClT 



L2 1 



C\\£, 



F22 = 2jTbf,{r~'bK,£tr-'S£t) < ^^^^^^JT{\\r-'dWEUh 
F23 = -2j6„(Et5wtT-i54^) +2Jr6;,(r-l5Etw^^r-l5£^^) 



5Vw^||^.||VE^||^. 



ClT 

IZ ClT IZ 

i^24 = -2J&„(5S^w^T-l8£,t) - 2j&„(S^-i,5w;;,r-i5£;J 

12 £2 t2 

i^25 = -2jb45ul£j;,T-'5£^) - 2j6,«-i,8£,t,r-i55^- 

,fe||2 IIV7!7fc||2 



< '-^\\T-^8V£t\\h 



2^2 



MfT 

ClT 



T-'5V£j:\\l. + %\5V^t\\Uys'-'\\h 

ClT 



/.iil2^^I|5Vu^IIl2||V4^L. 



/lllL2 



ClT 



By virtue of (2.1), F3 can be estimated as follows: 



71 A; 1 1 2 



< -i^„r||r-^5VE;;iiL. 



56k 



-(ll§^ 



2„,A;||2 



|5R^-i||2, 



ClT 



6 rci 



j/£||2 



|5E|;||L] 



The last two terms and i^s require no further manipulation and result in 



14 l/^T LZ 



fell 2 
L2 



Cl 



Il2 ; 



i^5<^||r-l5VE^||2, + l^||57^'^||2 



14C2 



ullL2 



12 " 



2o2 



/illL2 



ClT 



|8<I1l2 



Collecting all the estimates for ||t ^5VE)J||l2 and ||r ^5V£'^||l2, and using assumption (4.15), we get 



2C||E^ 



hllL2 



cj\\£j:\ 



L2 



T < V^T 



'2Cj\\£, 



hllL2 



r < ClT. 



/ji/2 ;ji/2 - ^1/2 

These conditions allow for cancellation of the problematic terms Fi^ and F22 with the fifth and sixth terms on the 
left hand side of (4.14). Finally, summation of the energy identity (4.14) and application of Gronwall inequality 
lead to ( [4T6| . □ 



Remark 4.2 (Smallness assumption). The error estimate (4.6) shows that (4.15) is actually a condition of the form 
(4.17) <c,<oo 

for a small enough constant Cs- This requirement is not a special characteristic of our method but rather a 
recurrent feature in the analysis of schemes for the Navier-Stokes equations. See, for instance [T^ I14| . 



4.3. Error Estimates for the Pressure. The control on the derivatives of the velocities provided by Theorem |4.2| 
enables us to obtain error estimates for the pressure. To do so, it is crucial that the discrete spaces are compatible 
in the sense of (2.11 ). This is the idea behind the following result. 



Theorem 4.3 (Error estimate for the pressure). // (4.2) and (4.16) are valid, then the following estimate holds 
(4.18) 



||e11,2(i2) <c(t + /i') 
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Proof. As already mentioned, the approximation errors and are actually solutions to (3.2 1 with a special 
right hand side composed of consistency terms. Condition (2.11) then allows us to write 



right hand 
(4.19) 
where 



l|e^lU^(0)< sup ^^<tB., 



Bi 
B3 



B2 = v„ sup — 1^ 1— 



sup 7 7 

sup 1^ 1^ , i34 = 2Vj. sup 

B5 = sup 7 7 , 

suffices to provide suitable bounds for each one of these terms. 

rlilv Vifi-\rp fm^ Ri finrl that 

ifc 



Bq ~ sup 



S2<||Ve;;||l.. 



So that it suffices to provide suitable boui 
We readily have, for Bi and B2, that 

i(i — sup 1^ 1^ ^ r 0±L^ 2 , 

Identity (4.9 1 can be used to express the numerator of B3 as 

B3 < sup II 1^ h sup n n h sup —5 



5 



sup 



n — n 1- sup n — n = yBsi. 

Inequality (2.6) and the regularity assumptions (4.1 1 imply 

6/i(5u'=,u'=,'L>/i) J, 
£31= sup r — 7 < ||u''||h2||ou''||l2 < ||5u'^||l2. 

the stability (2.13) of the projectors and the regularity 



To bound -632, -B33 and B35 we use inequality (2.15), 
assumptions (4.1), 

b,iu'^-\ S\v,) ^ , ^ , 



B32 = sup 



B33 = sup 
B35 = sup 



bh{E'^-\ulvh) 



< I|VE^-1||l2 



The first inequality in (2.17) yields 



^hjUf^ ,E,V^h) ^ II fc_i|||| fell 



<I|ve;;||l2. 



-B34 = sup 

conclusion, we have proved the bound 

li^sl < I1§u"||l2 + ||VS^-||l2 + ||VS^-1||l2 + HVEJ 



In 



ij—n ||E^ IIl2||VE,J1l2. 



li^sl <I1§u'=||l2 + ||VS^-||l2 + ||VS^ niL' 
Integrating by parts as in (2.1), we infer that 



;^1||l2+/.-V^||E^1||l2||VE;;||l2 



l|VEfj|L2. 



B^ = 2vr sup ^ \ 1^ sup ^ 7 7 ^ < ||8w''||l2 + ll'^^ ||l2 

"hSVh W^hWul veHi(n) II'^'iIIhi 



Finally, we see that 

= sup 



""f?''^ <r-^||5Sl|,2, B,^ sup „ „ 



;fc 

'^'-'F II Tl 



?^<IKIIl2. 
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It suffices now to realize that all the bounds involve consistency, interpolation or approximation errors and that 
they all have the right order. This concludes the proof. □ 

5. A Second Order Scheme 



Let us present a second order scheme for the solution of (1.51 and show its stability properties. We work in 
the setting of [2.2 and ^2.3 We first recall a three-term recursion inequality originally shown in [T3], which is 
instrumental to show stability. 

Proposition 5.1 (Three term recursion). The three term recursion equation 
(5.1) 3x''+^ - Ax'' + x''-'^ = yk > 1, 

has the following general solution 

1=2 s=2 

Let {y'^}fe>o be the solution to the three term recursion inequality 

3/'+i-4/ + /^i </+\ Vfc>l, 



with initial data y^ and y^ . If {x'^Jkyo is the solution to (5.1) with initial data = y^^ and x^ = y^ , then the 
following estimate holds 

y" < x\ \lv > 0. 
For {?/'^}j,>Q let 5^2/'^ denote the second order backward difference, i.e. 

52 / = ^(3/+i - 4/ + Vfc > 2. 

Let us now describe the scheme. We begin with an initialization step, in which we set 

(UtP^W^) = (utptw^), fc = 0,l. 

In other words, we compute the Stokes and elliptic-like projections of the initial data and the solution on the 
first time step. This initialization is only for ease of presentation as it clearly requires knowledge of the exact 
solution. In practice one can compute the projection of the initial data and then perform one step with the first 
order scheme of Section [3l 

We march in time, ior k ~ 2, . . . , K , as follows: 
Linear Momentum: Find (U^, Pjl) e V,, X Qh that solves 

(5.2a) + j.^(VUtVw^) + bh{Vl-\XJlvh) - [Pt.divVh) = 2iyr{cmlWf'\vh) + {f^Vh), 

(5.2b) {qh,divUt)=0, 

for all Vh G V/i, qh G Qh, where, for a time-discrete function 0"^, we introduced the second order extrapolation 

(5.3) c/)*^'* = 2/-1 
Angular Momentum: Compute £ Yh, solution of 

(5.4) jf ^^,a;J + c,{VWlVu:,,) + jb,,{VlWlu;h) + 

+ C2{divW''^,divuJh) + ii'riW^h^'^h) = 2iy^(curiUtu;,,) + (g^a;h) ■ 

for all uJh G V/i. 

This scheme turns out to be almost unconditionally stable, as shown in the following result. To avoid irrelevant 
technicalities, we assume that = g"^ = 0. 

Theorem 5.1 (Stability of second order scheme). Assume that the time step satisfies 

(5.5) r < ^. 



Then, the sequence {U^, W^, C [V/i]^ x Qh, solution of (5.2a|~(5.4), satisfies 
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where the constant depends on the material parameters and the values o/ {U^, W^, P^} for k = 0,1, but does not 
depend on the discretization parameters. 



Proof. We combine the techniques used to prove Proposition 3.1 and Theorem 5.1 of jT3]. We begin by setting 
Vh = 4rU{^ in (5.2a) and ujh — 4tWJJ in (5.4 1 and adding the resuh. Using (5.2b), we obtain 



3/ -4/-1 ' "■'=-2 



r - + 25(||5u^||^, +j||5w^;ii^.) + Wb'vUh +jI18'W^IIl= +4r K||VU^||^, + ci||VW,t||^.) 

+ 4c2r||divW;;||2, +16;.,r||W^||2, =8^..r(curiUt2W;;-52w;;), 



where 



2/^ = l|u;;ll^.+j|lw;;|l^. 



fc|i2 



Here we used the identity 



2a'=^(3a'= - 4a'^'-i + a''-^) = 310*^12 - 4|a'=-i|2 + |a'^-2|2 _^ 25|5c 



fe|2 



and, to produce the right hand side, we integrated by parts using (2.1 1 and employed the equahty 



20*^ 



;2 /fc 



which is a consequence of (5.3). Using (2.2 1 we obtain 



8i.,T(cur7Ut2W,';-52w;;) < 16i.,r|| W,';||^. + 4^.,,r|| VUf;!!^, + jIIS^W;;!!^ , 



IfH^i 
J 



IVU 



fc||2 



Since — v + assumption (5.5) yields 



2^2 



J 



= 4i/r 1 - 



> 2VT . 



The estimates of Proposition |5.1| imply the assertion 



□ 



Remark 5.1 (Time step constraint). Notice that the constraint on the time step (5.5 1, necessary for stability, is 
meaningful. First of all, the quantity on the right hand side has units of time. In addition, it is consistent with the 
fact that, for the classical Navier-Stokes equations (that is Vr = 0) no constraints are necessary for the stability of 
a second order semi-implicit discretization. 



6. Numerical Validation 

We now present a numerical validation of our error estimates. The implementation has been carried out with 
the help of the deal. II library, see We use the lowest order Taylor-Hood elements, that is Q2/Q1J so that 

1 = 2. The arising linear systems have been solved with the direct solver UMFPACK®. 

Consider a square domain Vl ~ (0, 1)^ C M^, and a smooth divergence-free linear velocity, pressure, and angular 
velocity defined by 

u(a;, y, t) = (sin(27ra; -|- t) sin(27ry -|- t), cos(27ra; -|- t) cos(27r2/ -|- t)y , 
P{x, y, t) = sin(27r(a; - y) + t), 
w{x, y, t) = sin(27rx -I- t) sin(27ry + t). 

To verify the £^(H^(i7)) error for the velocity and the £^(L^(J7)) error for the pressure we fix the relationship 
T = h"^, and consider a sequence of meshes with h = 2~* for 2 < i < 6. The corresponding errors are displayed in 
Figure |4j thereby showing clearly the predicted convergence rates. 

To validate the €°°(L^(f2)) error of the velocities we fix the relationship t = h^, and consider the same sequence 
of meshes. The corresponding errors are depicted in Figure [5] and exhibit the expected optimal rates. 
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Figure 4. i"^ (H.^ error of the velocities and £'^{L'^{il)) error of the pressure with respect to 
mesh size. The axes are in logarithmic scale. 




h 

Figure 5. £°° (L'^ (il)) error of the velocities with respect to mesh size. The axes are in logarithmic scale. 

7. Conclusions and Perspectives 

We have presented a first order, fully discrete semi-implicit scheme for the MNSE which is unconditionally stable 
and possesses optimal convergence rates in time and space. The scheme is semi-implicit, therefore it only involves, 
at every time-step, the solution of linear systems. In addition, the equations of linear and angular momentum 
are decoupled, which makes the implementation simpler and the scheme more efficient. To further decouple the 
unknowns, fractional time-stepping techniques can be incorporated, and we believe that their analysis shall not 
present difficulties beyond those already encountered in this work. 

We have also presented a formally second order scheme which is almost unconditionally stable and shares similar 
properties to the first order scheme, i.e., it is semi-implicit, decouples the linear and angular velocities and it can 
be easily simplified further with fractional time stepping techniques. The error analysis of such a scheme will be 
reported elsewhere, where in addition we will explore whether the stability condition is indeed a requirement of 
our scheme, or an artifact of our methods of proof. 
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The idea of pumping micropolar fluid through excitation of the spin equation was explored by testing a simple 
family of forcing terms g. It was observed computationally that the regimes of effective pumping and reverse 
pumping regimes are not well separated. In other words, very similar forcing terms g can induce very different 
effects in the velocity profile, or even opposite effects (reverse direction of the net flow). 

The most challenging extension of this work is towards the solution of the equations of ferrohydrodynamics: 
the MNSE with (1.8) coupled with the magnetostatic equations. The design, analysis and implementation of a 
scheme for this problem requires techniques and ideas well beyond those presented here, but will allow for more 
interesting and realistic simulations. This is part of future developments. 
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